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ABSTRACT 


A new method for enhancing convergence rate of iterative algorithms for the numerical integration of 
systems of partial differential equations has been developed. It is termed the Distributed Minimal 
Residual (DMR) method and it is based on general Krylov subspace methods. The DMR method 
differs from the Krylov subspace methods by the fact that the iterative acceleration factors are different 
from equation to equation in the system. At the same time, the DMR method can be viewed as an 
incomplete Newton iteration method. The DMR method has been applied to Euler equations of 
gasdynamics and incompressible Navier-Stokes equations. All numerical test cases were obtained 
using either explicit four stage Runge-Kutta or Euler implicit time integration. The formulation for the 
DMR method is general in nature and can be applied to explicit and implicit iterative algorithms for 
arbitrary systems of partial differential equations. 

INTRODUCTION 

After linearization caused by the discretization, the systems of governing equations associated with , 
say, fluid flows are recast into the following linear system of algebraic equations 

AX = b (1) 

where X is the vector of unknowns and A is an NxN matrix which depends on the discretized scheme, 
and is assumed to be non-singluar. The matrix A is usually sparse and as N becomes larger, it is not 
economical to solve the system of equations directly. Instead, iterative methods are usually utilized. 

The Conjugate Gradient (CG) method and the Conjugate Residual (CR) method, are widely used for 
approximating the solution of the system (Huynh, ref. 1; Faddeev and Faddeeva, ref. 2). Both 
methods give the exact solution in at most N steps in the absence of round-off errors. However, the 



CG method and the CR method require the matrix A to be symmetric, positive definite. A large number 
of generalizations of these methods applicable to systems with a non-symmetric matrix have been made. 
The success of the generalization of the CG and CR methods is reflected in the introduction of a series 
of algorithms capable of treating non-symmetric problems (OTHMIN by Vinsome, ref. 3; ORTHDIR 
and ORTHRES by Young and Jea, ref. 4; GMRES by Saad and Schultz, ref. 5; Wigton et al., ref. 6). 
The Minimal Residual method (Hafez, ref. 7) and the Generalized Nonlinear Minimal Residual method 
(Huang and Dulikravich, ref. 8) can be thought of as generalizations of the conjugate residual method. 

In this paper, a new method of enhancing convergence rate of iterative algorithms for systems of partial 
differential equations is developed. The method is entitled Distributed Minimal Residual (DMR) 
method (Lee et al., ref. 9-14) and it is related to a general Krylov subspace method from which it 
differs in two aspects. First, the DMR method attempts to improve on a straight application of a Krylov 
subspace method by using a separate sequence of acceleration factors for each equation in the system. 
In application of the DMR method to Euler equations of inviscid gasdynamics, for example, the 
acceleration factors for continuity equation differ from those for two momentum equations and for 
energy equation. This approach requires fewer consecutive solutions to be stored. Effectively, the 
DMR method periodically preconditions the system. Second, the DMR method does not involve the 
orthogonalization procedure which most of Krylov subspace methods utilize to reduce the number of 
numerical operations. The DMR method uses corrections from only two or three consecutive solutions 
for a successful application. 

The prime objective of this paper is to develop the theory of the DMR method and to examine the 
effectiveness of the DMR method by applying it to different systems of partial differential equations: 
Euler equations of inviscid gasdynamics and incompressible flow Navier-Stokes equations. Runge- 
Kutta time stepping method and Euler implicit method were used as two basic iterative algorithms. 

DISTRIBUTED MINIMAL RESIDUAL (DMR) METHOD 


Let us consider a system of partial differential equations that are integrated iteratively so that their 
residual vector at iteration level t is given by 


w t_9E^ dT^ 

dx + dy + dz 


( 2 ) 


where E l , F l , G l are the generalized flux vectors (at iteration level t) that act in the directions x, y, z, 
respectively. The future residual at iteration level t+1 is given by 
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( 3 ) 


D t + i 3E t+1 3F t+1 3G t+1 

R = 3 — + 3 - + 3 — 

ox dy dz 

Assume that each component of the solution vector at iteration level t+1 is extrapolated from the 
corresponding previous M consecutive iteration levels. Then, we can say that 

t+l 12 * 2 M A M 

q 1 = qj + co 1 Aq 1 + co^qj + ••• + co j Aqj 

q 2 = q 2 ^ 2 A ^2 ®2 A< ^2 

q l L =% + °i A % + °t A %+ - + m l Ac 1 l (4) 


Here, the subscripts 1, 2, 3, ..., L designate the particular component of the solution vector Q, that is, 
the particular equation in the system. The superscripts 1, 2, 3, ..., M designate the particular iteration 
level counting backward from the present iteration level, t. Thus, the superscript 1 means the first 
previous iteration level. The superscript 2 means the second previous iteration level, etc. This can be 
expressed in a more compact form as 
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( 6 ) 


Here, co's are the acceleration (weighting) factors to be calculated, A's are the iterative corrections 
computed with the original non-accelerated scheme, M denotes the total number of consecutive time 
steps from which the corrections are combined. 
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Using Taylor series expansion in time for R t+1 and truncating the terms that are higher than second 
order in At, Eq. 3 becomes approximately 


R t+1 = R 


M 

m=l 
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dy 
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3z 
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(7) 


The global domain residual can be defined as 

R l = ^ R tT R C (8) 

D 

where £ denotes summation over the computational domain D, and the superscript T represents 
D 

transpose of a vector. In order to minimize the future global residual, R t+1 , the CD's are determined 
from the following conditions 


3R t+1 
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= 0 


(9) 


From Eq. 8 this leads to 
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where 
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and 8 pr is the Kronecker delta. However, from Eq. 6 we have that 
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( 12 ) 
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Noticing that is not a function of co, it follows that 
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Then Eq. 13 becomes 
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For simplicity, let 
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Then, the system of algebraic equations (Eq. 15) can be written as 


M L 


EE 

n q 


(13) 


(14) 


(15) 


(16) 

(17) 


(18) 


263 



or 


r l 1 
C 11 

41 

... cjj 

-21 

C 11 

... Cli 1 

-11 

c 12 

4i 

••• 4.2 

r 21 

c 12 

-Ml 
• • • M2 

pit 

C 1L 

4l 

r lt 
... Cll 

-21 

C 1L 

••• 4L 1 

-12 

C 11 

-12 

C 21 

-12 
— C L1 

r 22 

C 11 

... cL{ 

-1M 

C 1L 

,1M 

C 2L 

,1M 
... Cjx 

-11 

C 1L 

4l m 



®2 


®L 

®1 


c ^ 1 


b 1 ! 

4 

4 

b? 


(19) 


representing the system of LxM linear algebraic equations for the LxM optimum acceleration factors co. 
For example, if we are periodically to combine corrections from M = 2 consecutive time steps to 
extrapolate the solution and to solve a system of L = 4 partial differential equations, we need to solve 
simultaneously LxM = 8 algebraic equations for 8 values of co. 

Notice that when the convergence is achieved, the b's become zero (Eq. 17), thus making the co's zero. 
In other words, the accuracy of the fully converged solution will not be affected by using the DMR 
method. Furthermore, if the matrix c™ is positive definite, it can be shown easily that the co's 

minimize the global residual, R l+1 , at iteration level t+1. Using a different sequence of acceleration 
factors for each partial differential equation in the original system is equivalent to using a different time 
step for each equation or selectively preconditioning the system. The DMR method, therefore, can be 
understood as the combination of a preconditioning method and a Krylov subspace method. Also, we 
can think of the DMR method as an incomplete Newton iteration. This point can be illustrated by the 
following fact. When the acceleration factors vary not only from equation to equation, but also from 
grid point to grid point, and when we use just one solution in the DMR formulation, it can be shown 
that the DMR method is equivalent to the Newton iterations. 


APPLICATION OF THE DMR METHOD 
TO EULER EQUATIONS OF GASDYNAMICS 


The introduction of the successful numerical algorithms such as the Euler implicit method and the 
explicit Runge-Kutta time stepping method made it relatively inexpensive to perform the numerical 
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integration of the systems of partial differential equations governing compressible flows. Most of such 
algorithms, however, suffer from slow convergence at low Mach numbers. The reasons for this are the 
rapidly increased stiffness and the singular' behavior of the original system of compressible flow 
equations at low Mach numbers. The singular behavior of the system near Mach number zero can be 
removed by eliminating the singularity of the system by a perturbation technique (Briley, ref. 15; Choi, 
ref. 16). The stiffness of the system at low Mach numbers can be reduced by preconditioning the 
system (Turkel, ref. 17; Choi, ref. 16). The DMR method is used to alleviate the difficulty associated 
with the increased stiffness of the Euler equations for low Mach number compressible flows. 


Euler Equations for Compressible flows 


The Euler equations for a two-dimensional unsteady inviscid flow expressed in a generalized non- 
orthogonal curvilinear coordinates (£,, q) without body forces or heat transfer, can be written in a vector 

form as 


= 0 

dt dt, 3q 
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where 
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The subscripts x and y represent first (partial) derivatives with respect to x and y, respectively. Here, p 
is the density, p is the thermodynamic pressure, e is the total energy per unit volume, while u, and v are 
the Cartesian velocity components along x and y axis, respectively. J is the Jacobian determinant, 

d(t'H) 


3(x,y) 


while U and V are the contravariant velocity vector components defined as 
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Numerical Algorithm 


The artificial dissipation suggested by Steger and Kutler (ref. 18) was used in the form 


D ( J Q) = -r- v4 [jQ] (23) 

UT A * ~ a 


where V 4 is the biharmonic differential operator in r\ coordinates and e is a parameter. The residual 

A 

vector R of Euler equations for compressible flow including the artificial dissipation is 


A 

R = 


3E | 3F | e 
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(24) 


After discretization, the governing equations become a set of ordinary differential equations, which can 
be integrated by the Runge-Kutta time stepping method (Jameson et al., ref. 18). 

Q° = Q l a 

AQ k = - a k AtR k " 1 k = 1,2,...,K (25) 

Qt+i = Qt + AQK 

where a k are the coefficients for each of the K stages of the Runge-Kutta scheme required to advance 
the solution from the time level t to the time level t+1. For example, a k = 1/4, 1/3, 1/2 and 1 for the 
four stage Runge-Kutta scheme. 


The time steps for each direction are estimated (MacCormack and Baldwin, ref. 19) from 
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(26) 


where c is the local speed of sound and CFL is the Courant-Friedrichs-Lewy number. The maximum 
time step is given as 


At^Atg 

AI “ At^ + At^ 


(27) 
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The implicit characteristic boundary procedure of Chakravarthy (ref. 20) was used, though the scheme 
itself is explicit. Entropy per unit mass (s = p/p 7 ), total enthalpy per unit mass, h = (e+p)/p, and flow 
angle (tan(a) = v/u) are specified at the inflow boundary. For a subsonic downstream outflow 
boundary (£, = constant), the equation corresponding to the negative eigenvalue, U - c(^ x + % ) , is 

substituted with a constant back pressure, p b . For a solid wall boundary (ri = constant), the equation 

2 2 1/2 

corresponding to the positive eigenvalue, V + c(\+ rj y ) , is substituted with a tangency boundary 
condition, V = 0. 


Upon applying the DMR method to the system of Euler equation of gasdynamics, Eq. 14, becomes 
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where A and B are the Jacobian matrices in the transformed coordinates 
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dE 

dQ 



(29) 


Results for Compressible Euler Equations 

A two-dimensional flow analysis code has been developed in FORTRAN according to the previous 
theory for Euler equations of gasdynamics using finite differencing. All computational results were 
obtained on CRAY-YMP at NAS facility using automatic vectorization. 

The test case for the code was flow around a circular cylinder. The outer boundary of the computational 
domain was located at 20 times the radius of the cylinder. A 66x32 cell computational grid was used in 
this test case. The computations were performed with and without the DMR method in conjunction with 
the four stage Runge-Kutta (RK) scheme. The convergence histories are plotted in terms of the number 
of iterations and in terms of the CPU time (Fig. 1). The maximum allowable CFL number (CFL = 2.8) 
was used in both accelerated and non-accelerated computations. The free stream Mach number was 
chosen to be 0.05 which is practically an incompressible flow. The DMR method saves over 60 % of 
total CPU time in this critical flow test case. The surface pressure coefficient (Fig. 2) matches well with 
the incompressible potential flow solution. 
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Figure 1 Convergence histories for the inviscid flow around a circular cylinder with = 0.05 

Wall Pressure Coefficient 



Figure 2 Wall pressure coefficient distribution for the inviscid flow around a circular cylinder: 
numerical with = 0.05 (solid line); analytical with = 0 (dotted line) 


APPLICATION OF THE DMR METHOD 
TO INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 

The main difficulty associated with the incompressible flow computations is caused by the absence 
of a time derivative term in the continuity equation. One of the methods for solving the incompressible 
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Navier- Stokes equations was originated by Chorin (ref. 21). In this concept, an artificially time 

dtn/Rl 

dependent derivative term - gp- is added to the continuity equation with a user specified control 

parameter p. The artificial time derivative diminishes as the solution converges to its steady state. The 

added term forces the system to be of a mixed parabolic -hyperbolic type, which allows the use of time 
marching techniques. Later, Choi and Merkle (ref. 22), and Kwak et al. (ref. 23) used an Alternating 
Direction Implicit (ADI) method in conjunction with the artificial compressibility method. 

Incompressible Flow Navier-Stokes Equations 

The two-dimensional Navier-Stokes equations in a general non-orthogonal curvilinear coordinates r) 
are given as 


f + f + f = D2(JQ) 

drj 


The solution vector and the flux vectors in the transformed coordinates are given as 
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where p is the pressure. Notice that the artificial compressibility has been added in the continuity 
equation. The physical viscous terms in the general coordinates are given by 

£> 2 (JQ) = [jgjj(jQ) j] . (32) 

where gy is the contravariant matrix tensor 

gij = Vx'iVx'j (33) 

Here, x) means ^ or q depending on the index i 

S =“^diag(0,l,l) (34) 

where Re is the Reynolds number. 
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The Navier-Stokes equations are mixed parabolic/hyperbolic partial differential equations. According to 
the eigenvalue analysis of the hyperbolic part of the equations, the Jacobian matrices in the transformed 
coordinates have real eigenvalues 


A=^ = K(u,U; y ) 
BQ 

B = — = K(V ,rq x ,ri )) 

BQ 


where the matrix K is defined as 
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(35) 

(36) 


(37) 


Here, kj and k 2 are either and or r| x and T) y depending on the direction to be considered, and k = 
kju + k 2 v. The eigenvalues of the matrix K are given by 


A = diag(k - c, k + c, k) 


(38) 


where the equivalent speed of sound, c, is given as 


c = 



k 2 + |3(kj + k\) 


(39) 


Notice that one of the eigenvalues is negative. This means that the incompressible flow is equivalently 
"subsonic" in the sense of different signs of the eigenvalues and that c will influence stiffness of the 
system. Thus, the direction of characteristics should be considered when applying boundary 
conditions. 
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Numerical Methods 


The residual vector including the fourth order artificial dissipation (Eq. 23) is defined as 


A r)F dF 1 ~ ~ 

R ^ + r^ 2 ( J Q) + D (JQ) 

d^ dp 


(40) 


After spatial derivative terms were discretized, the governing equations were integrated either by the 
explicit Runge-Kutta time-stepping algorithm (Eq. 25) or by an Euler implicit method with approximate 
factorization (Beam and Warming, ref. 24). To reduce the computational effort, the artificial dissipation 
and the viscous part of the residual are calculated only once every global time level and kept unchanged 
during the four stages of the Runge-Kutta scheme. This does not deteriorate the stability of the time 
stepping algorithm. 


The Euler implicit scheme with factorization for the incompressible Navier-Stokes equations results in 
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Time Step Limitations and Boundary Conditions 


The allowable time increments of the explicit scheme are severely restricted by the stability limit, while 
for an implicit scheme the time step restrictions are caused by the factorization errors. The time step is 
determined by considering the hyperbolic part of the system and the parabolic part of the system 
separately and by combining these time steps as suggested by MacCormack and Baldwin (ref. 19). The 
system becomes hyperbolic when viscosity is neglected. Then, the stability bound of the resulting 
system is determined by the CFL (Courant-Friedrichs-Lewy) number. The maximum allowable time 
steps for each of the coordinate directions are defined as 
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so that the combined maximum time step for the hyperbolic part of the system is defined by 
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When the convective part of the acceleration is neglected, the system becomes of parabolic type. The 

stability of the parabolic type system is dictated by the non-dimensional number a (von Neumann 
number). For each generalized coordinate direction, the maximum time steps are defined by 


* aRe 

A te = ITT 


, aRe 
Atp11_ feT 


(44) 


and the combined maximum time step for the parabolic part is given by 

At At P^ A Vl 

^ At p ^ + At pT) 

The total maximum time step is estimated conservatively as 

Al= At - At P 

At h + At p 


(45) 


(46) 


For the explicit Runge-Kutta method, Eq. 46 was used to estimate the maximum time step. However, 
for the Euler implicit method, only CFL limitation was used to compute the time step, that is 


At = 


At)^At hr[ 
At h ^ + At hT1 


(47) 


It was assumed that the flow is inviscid at the inlet and exit planes causing the system of equations to 
become hyperbolic in time near the inlet and exit. As stated earlier, the incompressible Navier-Stokes 
equations have one negative eigenvalue, and the rest of the eigenvalues are positive. Thus, one 
equation should be considered with two boundary conditions at the inlet. At the exit, two equations 
with one boundary condition must be applied. At the inlet, u and v velocity vector components were 
specified, while the back pressure p was specified at the exit. The flow was assumed to be locally one- 
dimensional at the inlet and exit boundaries in order to transform locally the equation into the 
characteristic form. At the solid wall, the velocity vector components, u and v, were set to zero, and 
the surface pressure was extrapolated from the grid points next to the wall from the condition that^ = 



Residual Smoothing 


One of the successful attempts to accelerate the convergence of the Runge-Kutta scheme is Implicit 
Residual Smoothing (IRS) introduced by Jameson and Baker (ref. 25). With this method, it is possible 
to use much higher values of CFL. The residual is smoothed through the following equation 


i - esJTi 
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where 8 2 designates the central difference operator for a respective second derivative, and 0 is the 
smoothing coefficient. Thus, when using the IRS we have to solve two scalar tri-diagonal matrices. 
Since their coefficients are constants, the tri-diagonal matrices are decomposed into upper and lower bi- 
diagonal matrices so that at every application of the IRS only forward and backward substitutions are 
needed to get the smoothed residual. 


The application of the DMR method to incompressible Navier-Stokes equations differs from the 
formulation for its application to the Euler equations of gasdynamics only by the following term 
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Computational Results for Navier-Stokes Equations 

A steady, laminar, viscous flow normal to a solid wall (Hiemenz flow) was the first test case. Reason 
for this choice of the test case is that the analytic solution for the Hiemenz flow is known (Panton, ref. 
26). The accuracy of the codes (the explicit Runge-Kutta method and the Euler implicit method) can be 
verified by comparing the computed solution with the analytic solution. 

The flow corresponding to the Reynolds number 400 based on the free stream velocity and a body 
dimension, R 0 , of the wall was computed with and without the DMR method in conjunction with 

explicit and implicit codes. The computational grid consisted of 60x29 cells, and the dimensions of the 
computational domain were H = R 0 and L = 2R 0 . In the case of an explicit Runge-Kutta (RK) method, 
the maximum allowable CFL number of 2.8 was used and the von Neumann number was a = 0.4. A 
small amount of the fourth order artificial dissipation was added to get a smooth solution (e = 0.05). 
Using numerical experimentation it was found that the fastest convergence is obtained with the artificial 
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compressibility coefficient (3 = 2, and that the DMR method should be applied every 10 iterations by 
combining 3 consecutive solutions. 



Figure 3 Distributions of wall surface velocity gradient for Hiemenz flow 
(RK: solid line; analytic solution: circles) 



The computed distribution of the wall surface velocity gradient, was compared with that of the 
analytic solution (Fig. 3), showing an excellent agreement. Figure 4 shows that the residual was 
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reduced 12 orders of magnitude in 5000 iterations without the DMR method, while the same reduction 
in residual could be achieved in 2000 iterations with the DMR method indicating 60 % reduction in 
CPU time. The implicit residual smoothing was also implemented with and without the DMR method. 
The basic RK method gives the slowest convergence, the IRS gives faster convergence than the basic 
RK method, while the DMR method gave the second best convergence. The most rapid convergence in 
terms of the number of iterations was achieved by combining the implicit residual smoothing and the 
DMR method. However, the DMR method alone offered maximum time savings (over 55%). 



Figure 5 Distributions of wall surface velocity gradient for Hiemenz flow 

(Euler implicit: solid line; analytic solution: circles) 



Number of Iterations Cpj Time (sec) 


Figure 6 Convergence histories of the Euler implicit method for Hiemenz flow with Re = 400 
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The implicit code was also exercised for the test of the Hiemenz flow with the same conditions as in the 
test case for the explicit code (R.e = 400). The computed surface velocity gradient distribution was 
compared with the analytic solution (Fig. 5). Good agreement can be observed. CFL number of 10 
was used in this computation. Also, the fourth order artificial dissipation with e = 0.25 was added. 
The optimal value of the artificial compressibility coefficient p was found by numerical experiments to 
be P = 5. The DMR method was found to give the fastest convergence when applied to the implicit 
Euler scheme every 5 iterations by combining 5 consecutive solutions. Figure 6 shows that the DMR 
method offers approximately 60% reduction in CPU time indicating that the DMR method can be 
successfully applied to implicit methods. 



Nunber of Iterations Cpu Time (sec) 

Figure 7 Convergence histories for viscous flow around a circular cylinder with Re = 20 (RK) 

The last test case was a laminar flow around a circular cylinder. The highly clustered grid of 66x44 
cells was used. Flow with Reynolds numbers of 20 was computed with the RK method and the Euler 
implicit method. The CFL and von Neumann numbers were CFL = 2.8 and a =0.4, respectively, for 
the RK method, and CFL =10 was used for the Euler implicit method. The DMR method was applied 
every 30 iterations for the RK method, while every 10 iterations for the Euler implicit method. For 
both methods, two consecutive solutions were used with the DMR method, though these combinations 
of the number of solutions and the frequency of the DMR application are not optimal. The artificial 
compressibility coefficient was P = 1 for both methods. The convergence histories of the RK method 
and the Euler implicit method with and without the DMR method are presented in Figures 7 and 8. The 
DMR method offers more savings with the Euler implicit method (30%) than with the RK method 
(10%). The wall pressure distributions and the wall vorticity distributions were compared with the 
computational results of Choi (1989) in Figures 9 and 10 showing reasonable agreement. 
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Convergence History 



Figure 8 Convergence histories for viscous flow around a circular cylinder with Re = 20 
(Euler implicit) 


Wall Pressure Coefficient 



(a) Pressure coefficient 


Figure 9 


Wall Vorticity 



Angle (deg.) 


(b) Vorticity 


Wall pressure coefficient distributions and vorticity distributions for flow 
around a circular cylinder at Re = 20 (RK) 
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(a) Pressure coefficient (b) Vorticity 


Figure 10 Wall pressure coefficient distributions and vorticity distributions for flow 
around a circular cylinder at Re = 40 (Euler implicit) 

CONCLUSIONS 

The DMR method was found capable of reducing the computation time by 20-80% depending on 
the test case. When directly compared with an implicit residual smoothing, the DMR method performed 
consistently better and more reliably. 
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